%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mapping reduced form coefficients into structural coefficients 
% Exact identification for skill bias project
% this version: October 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% inputs:
% BB                reduced form coefficients for ndraws
% Omega             VCV from reduced form residuals for ndraws
% nlags             number of lags
% ndraws            number of draws
% HORIZON_ident     identification horizon
% c                 choice of constant in the reduced form VAR
% num_vars          number of variables in the VAR
% nshocks           number of identified shocks
% hndx              index of hours in the VAR
% do_level          hours in levels of first differences

% output:
% Response          Responses to one standard deviation of all shocks for all horizons and all draws
% Atilde            matrix A from the structural representation
% Respres           Response to (all) one standard deviation residual shock(s)
%                   size is num_vars,HORIZON_ident,ndraws

function [Response,Atilde,Respres]=ident(BB,Omega,nlags,ndraws,HORIZON_ident,c,num_vars,nshocks,hndx,do_level)

Response = zeros(num_vars*num_vars,HORIZON_ident,ndraws);
Atilde = zeros(num_vars,num_vars,ndraws);
Respres = zeros(num_vars,HORIZON_ident,ndraws);

for n = 1:ndraws

    % Iteration
    dnum=num2str(n);
    if dnum(end)==num2str(0) && dnum(end-1)==num2str(0)
        disp('iteration'),disp(n)
    end
    
    % Step I: 
    % transform the reduced coefficients into structural coefficents Phi
    Btilde = zeros(num_vars,num_vars,HORIZON_ident);
    for k=1:nlags
        Btilde(:,:,k) = BB(:,((k-1)*num_vars+c+1):(k*num_vars+c),n);
    end;
    Phi = zeros(num_vars,num_vars,HORIZON_ident+1);
    Phi(:,:,1) = eye(num_vars,num_vars); % this is time zero, i.e. the impact response

    for s = 2 : HORIZON_ident,
        Phitemp = zeros(num_vars,num_vars);
        for j= 1:(s-1),
            Phitemp = Phitemp + Phi(:,:,(s-j))*Btilde(:,:,j);
        end;
        Phi(:,:,s) = Phitemp;
    end;

    % Step II: Implementing the restrictions
    % i.e. get CCtilde and A, the identification horizon is flexible here
    % and may be adjusted
    CCtilde = eye(num_vars,num_vars);
    for j= 2:HORIZON_ident+1
        CCtilde = CCtilde + Phi(:,:,j);
    end;

    % we are now looking for some matrix A, such that AA'=Omega and such that
    % CCtilde*A is lower triangular
    Sigma=CCtilde*Omega(:,:,n)*CCtilde';
    Qtrans = chol(Sigma);
    Q=Qtrans';
    Atilde(:,:,n) = CCtilde\Q;       

    % Step III: IRFs for all shocks (separately)
    for shock_counter = 1:num_vars,       
        shock_vector = zeros(num_vars,1);
        shock_vector(shock_counter) = 1;
        CCneutilde = zeros(num_vars,num_vars);
        for time_counter = 1 : HORIZON_ident,
            CCneutilde = CCneutilde + Phi(:,:,time_counter);
            Resphelp1 = CCneutilde*Atilde(:,:,n)*shock_vector; % for differences
            Resphelp2 = Phi(:,:,time_counter)*Atilde(:,:,n)*shock_vector; % for levels
            Response((shock_counter-1)*num_vars+1:(shock_counter-1)*num_vars+num_vars,time_counter,n) = Resphelp1;
            if do_level==1
                Response((shock_counter-1)*num_vars+hndx,time_counter,n) = Resphelp2(hndx,:)*100;
            end
        end;
    end;

    % IRFs for all residual shock (bunched together)
    shock_vector = zeros(num_vars,1);
    shock_vector(nshocks+1:num_vars) = 1;
    CCneutilde = zeros(num_vars,num_vars);
    for time_counter = 1 : HORIZON_ident,
        CCneutilde = CCneutilde + Phi(:,:,time_counter);
        Resphelp1 = CCneutilde*Atilde(:,:,n)*shock_vector; % for differences
        Resphelp2 = Phi(:,:,time_counter)*Atilde(:,:,n)*shock_vector; % for levels
        Respres(:,time_counter,n) = Resphelp1;
        if do_level==1
            Respres(hndx,time_counter,n) = Resphelp2(hndx,:)*100;
        end
    end;
    
end; % draws loop